clearvars -except b

load tprob_data 


M = max(tprob_m);
pxest2 = 99999+zeros(M,max(tprob_t));


yvals = [1:1:19]';

for m = 1:M
    for y=1:length(yvals);
    index = find(tprob_m==m & tprob_t <= y);
    if ~isempty(index)
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index)]);
    pxest2(m,y) = temp(2);
    end
    end
end



T = 7; beta = .95;

gammaae = zeros(M,19);
for m = 1:M
    %m
    for y=1:length(yvals);
    if pxest2(m,y) < 99999
    gamma = 1;
    B = 1;
    for t=2:T
    temp = beta^(t-1);
    temp2 = (beta^(t-1))*( pxest2(m,y)^(t-1));
    B = B+temp;
    gamma = gamma + temp2;
    end
    gamma = gamma/B;
    gammaae(m,y) = gamma;          
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load est_res_s2 rentd  M con_m t
load datafromstata gamma 

statmwtp = rentd;

picmwtpae = 0*rentd;  
for m = 1:M
    for y=3:length(yvals);
        index = find(con_m==m & y==t);
        picmwtpae(index) = statmwtp(index)/gammaae(m,y);
    end  
end

for m = 1:M
    index2 = find(con_m==m & t >= 3);
    res2(m,:) =  mean(picmwtpae(index2));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rob_re = res2;
save rob_re rob_re





